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Abstract 

In a planetary or satellite system, idealized as n small bodies in initially coplanar, concentric orbits 
around a large central body, obeying Newtonian point-particle mechanics, resonant perturbations will 
cause dynamical evolution of the orbital radii except under highly specific mutual relationships, here 
derived analytically. In particular, the most stable situation is achieved (in this idealized model) only 
when each planetary orbit is roughly twice as far from the Sun as the preceding one, as observed 
empirically already by Titius (1766) and Bode (1778) . Simplifying the problem by reformulating it 
as a hierarchical sequence of [unrestricted] 3-body problems, in which gravitational interactions are 
ignored except between the central body and the body of interest and the next outwardly orbitally 
adjacent body, it is proved that the resonant perturbations from the outer body will destabilize the 
inner body [& conversely] unless its mean orbital radius is a unique and specific multiple (3 of that 
of the inner body. In this way a sequence of concentric orbits can each stabilize the next adiacent 
inner orbit, until only the outermost orbit remains; but it is already tied to the collection of inner 
orbits by conservation of total angular momentum and so the entire configuration becomes stabilized. 
Let /x = Myiax/M denote the ratio of the mass Mmax of the largest small body to that, M, of the 
large central body; in our Solar System, [i is less than 10~ 3 . Expanding (3 in a power series in fi, 
the lowest-order terms for distal multiplier (3 and phase shift cf) are found to start with the following 
universal constants [for m = 1, 2, 3] 

f3 = /3+ = l/y / (3/2)(2/ 3 ) - 1 = 1.794980 + , <j> = m(vr/2) + .. 

which agrees with the observed distal ratios between Jupiter and the Asteroids, and between Saturn 
&: Jupiter, and Uranus & Saturn. In view of Kepler's Law this implies that the stabilizing resonances 
must be of the type 2:1 or 5:3 or 3:1 (to the extent that only 3-body interactions and only integers not 
exceeding 5 are considered). Only analytic complexity has inhibited the present theory from being 
extended to 4-body and 5-body interactions, as a result of which the observed resonant structure 
of the Solar System [Ovenden & Roy; Molchanov] can be understood dynamically as the inevitable 
result of resonant perturbational evolution. 

Keywords: Planets and satellites: general-solar system: formation 
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Introduction 



The distribution of planets in the Solar system is not random and their mean distance, from the Sun, 
when numbered from the centre, form a rough geometric progression: 

» 1.75 (1) 

or 

r n = 0.4 + 0.15 x 2 n ,n = -oo, 1.. .8 (2) 

This last relation is known as Titius-Bode law, and it roughly describes the planetary semi-major axes 
in astronomical units, with Mercury assigned i = — oo , Venus i — 1, Earth i — 2, etc. Usually the asteroid 
belt is counted as i = 4. The law fits the planets Venus through Uranus quite well, and successfully 
predicted the existence and locations of Uranus and the asteroids. However, (i) the law breaks down badly 
for Neptune and Pluto; (ii) there is no reason why Mercury should have % = — oo , rather than % — 0, 
except that it fits better that way; (iii) the total mass of the asteroid belt is far smaller than the mass 
of any planet, so it is not clear that it should be counted as one. When this relationship was discovered 
by Titius of Wittenberg in 1766 and published by Bode six years later, it gave good agreement with the 
actual mean distances of the planets that were then known-Mercury (0.39 AU), Venus (0.72 AU), Earth 
(1.0 AU), Mars (1.52 AU), Jupiter (5.2 AU), and Saturn (9.55 AU). Uranus, discovered in 1781, has mean 
orbital distance 19.2 AU, which also agrees. The asteroid Ceres, discovered 1801, has mean orbital distance 
2.77 AU, which fills the apparent gap between Mars and Jupiter. However, Neptune, discovered 1846, has 
mean orbital distance 30.1 AU, and Pluto, discovered 1930, has mean orbital distance 39.5 AU; these are 
large discrepancies from the positions 38.8 AU and 77.2 AU, respectively, predicted by Bode's law. It is 
thought that the outer planets do not fit as well due to millions of tears of comet impacts and a looser 
gravitational hold. Also, Pluto is thought to have been created outside our solar system and the law would 
therefore not apply to it. 

For the next 200 years there has been a debate in astronomy as to whether Bode's Law is simply chance 
or if it underlines a yet unseen force in the creation of the solar system. A fairly comprehensive history of 
the law and attempts to explain it up to the year 1971 can be found in Nieto [43]. Some theories of the 
origin of the solar system have tried to explain the apparent regularity in the mean orbital distances of 
the planets, arguing that it could not arise by chance, but must be a manifestation of the laws of physics. 
Some astronomers hold that the deviation of Neptune and Pluto from their predicted positions signifies 
that they are no longer at their original positions in the solar system. However, since Bode's law is not 
a law in the usual scientific sense, i.e., it is not universal and invariant, it alone should not be taken as 
evidence for such a conclusion. Under the assumption that planets need a certain amount of space to form 
without competing for material, random numbers subjected to computer simulations of early solar systems 
have been done which show planet spacing that follow similar laws [44]. Despite this, many scientists over 
the centuries have come up with explanations for the law. One theory, presented by Stuart Weidenchilling 
and Donald Davis, claims that the planets ended up where they are due to the frictional drag, bringing 
the smaller planets closer to the sun and "gravitational perturbations from embryonic planets" Under the 
right conditions there would be favourable locations for planets to form [45] . Another recent theory is that 
a simple numerical sequence that "combines scale and rotational invariance, the points at which a physical 
parameter, such as density, reaches a maximum or a minimum will always follow a simple relationship". 
Using this insight, Francois Graner and Berengere Dubrulle have improved Bode's Law to give better 
results. They claim that to believe in Bode's Law is to believe in scale invariance in the early solar system 
[46]. Summarizing most modern arguments concerning the validity of Bode's law can be assigned to one 
of three broad classes: 
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1. Attempts to elucidate the physical processes leading to the Bode law. These are based on a variety 
of mechanisms, including dynamical instabilities in the protoplanetary disk ([47], [48], [49]), gravitational 
interactions between planetesimals [44], or long-term instabilities of the planetary orbits ([50], [51], [52]). 

2. Discussions that ignore physics but try to assess whether the success of the Bode law is statistically 
significant ([53], [54], [55]). Conclusions go from the Bode law being real rather than artifactual to the 
contrary. As observed by Hayes & Tremaine [56], these last conclusions are flawed because of some 
assumptions in the analyses. 

3. Discussions of other laws that may influence the spacing of the planets. Many of these involve resonances 
between the mean motions of the planets, such as [57] (but see [58]), [59], [60]. 

For what concerns point 3, the commensurability and resonance phenomena of the solar system motion 
structure - (including planets, asteroids, planet's satellites) - have been object of detailed discussions and 
experimental examination in the last years. Systematic observations and measurements have been carried 
out using all available means. A considerable bulk of empiric data for the solar system has been obtained 
from the "Voyager 2" mission. The understanding of the inevitable resonance character of evolving mature 
oscillation systems leads to series of interesting concepts about the resonance character of the Solar system 
motion dynamics. One of the most outstanding among them is Molchanov's hypothesis about the complete 
resonance character of the large planets in-orbit motion. A.M. Molchanov has noticed that the mean motion 
of the nine large planets are related approximately by nine linear homogeneous equations , 

k{ j) m + + 4 j) n 9 = 0,j = 1,....9 (3) 

with integer coefficients: 

k?, ,k? (4) 

The mean motions of the Jovian, Saturn and Uranian satellites are related by similar equations. If asteroids 
are considered as planets and excluding Pluto, it is established [61] that the planetary distances obey the 
following regularity: 

— = 1.75 ±0.20 (5) 

a-k 

i.e. the relation of the semi-major orbital axes of neighboring planets is almost constant. 

In the present paper, we perform a dynamical derivation of Bode's law. The basic strategy in the 
following is as follows: the famous "Poincare map" which is now often implemented by physicists using 
numerical integration [34] will be implemented analytically. In the present context, this means that an 
arbitrary line-segment will be defined transverse to a periodic generating orbit, and nearby orbits will be 
followed until this 1-dimensional "surface of section" is intersected a second time. This defines a continuous 
mapping of the section into itself, of which map the initial periodic solution is now a fixed point. Behavior 
of nearby orbits can now be studied by considering the simpler problem of iterating this map in a small 
neighborhood of the fixed point. It will be proved below that no fixed point exists (i.e. no periodic solution 
exists) unless at ji = the distal ratio defined in the Abstract is precisely (3 . 
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Problem Formulation 



Let x e E 2 (i = 1, 2, 3, — , ri) denote 2- vectors or elements of 2-dimensional real Euclidean 
space E which represent the positions of n point-particles of masses M. in orbit around a central 

body of mass M and position x°. Let v s E 2 denote their velocities x 1 = dx/dt with respect to 
time t, where := d/dt. Let x-y = (x,y) denote the scalar product between any two vectors x, y, 

1/2 

and let llxll := (x,x) denote the norm (length) of any vector. 

It is assumed that u, := MIM « 1. Now rescale the masses so that M = 1 and the smaller 

i i o o 

masses are Lt « 1. Let u, denote the largest of the small masses, and rewrite them as 

n 'max b 

\i. = u,u° = \i-£°-\i MAX , < u. < 1, ( 0< e° < 1 ), (6) 

so that |i is the perturbation parameter. All that will be PROVED here is for |i 'sufficiently 
small', though there are good reasons for believing that an analytical continuation (in the manner 
of Poincare) can be made all the way from [i = to |i = 1; this is because the Leray-Schauder 
Index of the 'generating solution' at \l = will be proved to be unity, and the fact that this 
integer-valued topological invariant, specifying in some sense the true multiplicity of actually 
existing solutions, is a continuous homotopy invariant and so can change its value discontinuously, 
i.e. terminate', ONLY at a bifurcation or singularity of the solution (i.e. a collision or an 
ejection), demonstrates that the homotopy on |i e [0,1] is legitimate unless there is an intermediate 
value of (i < 1 at which there is a collision between two bodies or one body is ejected to infinity; 
this makes rigorous Stromgren's Principle of Natural Termination, discovered empirically by 
numerical integration [1], in which a family of periodic solutions being studied by the variation of 
a parameter can cease to exist only at an ejection/collision event. 

Now introduce relative coordinates, in which each x 1 is replaced by (x 1 - x°), although for 
convenience the latter will be renamed as x; it is well known that then the system becomes, for the 
indices i = 1, 2, ■■■ , n, 

X = V, (7a) 
n 

v 1 = -G-(l + tytf/wxh 3 ) - ^^[[xVllA 3 ] - [(x 1 - x)/^ - x 1 )!! 3 ]] . (7b) 

j=i 
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where the absolute acceleration of the coordinate system's origin, which now coincides with the 
position of the Sun or central body, is the source of the first term in the summation (cf. Kurth 
[24], pp. 84-85), and where G denotes the Newton-Cavendish parameter or gravitational constant. 

Now, following Kurth, assume that each planet is affected significantly only by the Sun and the 
next further outward planet: 

V k = -G(l + ^)-(//ll/l| 3 ) + G-U^m _ _ _ ^/H^H^j, 

k = O'j), j = (i + 1); m = m(k); m(i) = j; m(j) = i; (i = 1, 2, ■■■ , n - 1). (7c) 

Again, following Kurth, I shall solve the Newtonian system (7a,c) by successive approximations, 
but in a modified manner. Note that if = 0, then (7c) can be solved by concentric circular 
solutions, with frequencies Cfl k given by Kepler's Third Law (1619) as 

(co/ = G-(l + |a k )/(p k ) 3 , p k := h/h. (8) 

For the purpose of successive approximations, replace (7c) by the equivalent ordinary 
differential equation (ODE) 



v k = -(co k )V + G[l k 



m km k3 mm3 

[(x - X )I\\X - X II ] - [x l\\X II ] 



+ [(co k ) 2 - Gil + |i k )/n/n 3 ]-/, (9) 



which makes clear the fact that the system is just a perturbation of the harmonic oscillator 
problem. 

We shall need the following lemma, which is an obvious modification of well known basic results 
in the theory of ODEs [25], [26], [27], [28]. 

Lemmas from ODE Theory 

DEFINITION . An m-vector function f: E m -> E m has a global Lipschitz constant K if there is a 
positive real number K > such that, for all z e E m , 

ii f(z 2 ) -f{z) ii < k- ii z - z l H. (10) 

THEOREM . Consider the m-vector ODE Initial Condition Problem (ICP) 

z = f(z), z(0) = z°, ( < t < T), (11) 

where f satisfies a global Lipschitz condition, and where the maximum time T is arbitrary but 
finite. Let A denote an arbitrary m x m constant real matrix, and define the vector function g by 
g(z) '■= f(z) - Az, which also is globally Lipschitzian, with Lipschitz constant now /k + llAll], 
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where I LA II denotes the Euclidean norm of A. Now construct the sequence of functions { z (t) } by 
successive solution of the forced linear ODE 

z k = Az k + g(z k -\t)), z(0) = z°, ( k = 2, 3, - ), (12) 

where the initial iterate z\t) is any arbitrary continuous function of t on [0, 7] such that its 

initial value z\0) = z°; or, equivalently, construct the sequence by 

t 

z\t) = exp{Ai)-z + | exp(A[t - x])-g(z k '\x)) dx, < t < 7, (13) 

o 

i.e. by integrating the right-hand-side to obtain the left-hand side, and then inserting the result 
into the right-hand-side and repeating the operation. The resulting sequence is guaranteed to 
converge uniformly on [0, 7] to the exact and unique solution of (11), no matter what is the 
arbitrary initial iterate z l (t), and no matter what is the arbitrary matrix A, and no matter how 
large is the [fixed] time-interval [0,7]. 

PROOF . Any standard work on ODE's will show how to use Lagrange's method of variation of 
constants to prove that the ODE ICP 

z = Az + g(z), z(Q) = z°, ( < t < T ), (14) 

is completely equivalent to the Volterra integral equation problem 

t 

z(t) = exp(At)-z° + J exp(A[t - x])-g(z(x)) dx, < t < 7, (15) 

o 

where 

oo 

exp(At) := ^ L A k -t k ik!. (16) 

k=0 

Therefore, if the sequence in (12)-(13) converges uniformly, then the result satisfies (15) and 
so is the solution of (11). Because every continuous function on a closed and bounded subset of a 
finite-dimensional Euclidean space has (and assumes) a finite maximum and minimum, there exists (on 
the closed, bounded subset [0,7] cz E 1 ) 

4> = MAxi II z\t) - z\t) II }. (17) 
te [o,T] 

By repeatedly applying the triangle inequality to (16) one finds readily that 

oo 

\\exp{At)\\ < J^\\A\\ k -t k k! = exp{\\A\H) < y : = exp(\\A\\-T). (18) 

k=0 
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Accordingly it is easy to prove by induction that (setting k = K + llAll) 

ii z k+ \t) - z\t) ii < (j) o -(Y-K-O k /^/, 



(19) 



LIM { II z\t) - z\t) ii } < [ II z k+ \t) - z\t) H < ^exp(y-K { T), 



(20) 



k=0 



for each t e [0,7], so that { z if) } is a Cauchy sequence (for each t) and its limit z(t) must 
exist. Consideration of the approximating sums E N and X N+1 to the series in (20) shows that the 
convergence is uniform. □ 

In the sequel, we shall frequently use the 2 x 2 identity matrix 1^ and its 'imaginary' skew- 
symmetric counterpart J ^ = -{J ^ defined as 



7 = 

2 



1 , o 

, 1 



J = 

2 



, 1 
-1 , o 



(jy = -i. 



For present purposes, the result that 



exp{JA) = 



cos(§), sin(§) 
-sin(§), cos(§) 



(21) 



(22) 



will be all-important. (To prove it, insert (21c) into (16).) A similar result which will be 

needed is that if the 4 x 4 matrix A is defined by (using MATLAB notation) A = [0, 1^, -Cfl 2 7 2> 0], 
then it is easy to prove by induction that 

A 2k = (-l) k -co 2k -7 , 



A 2k+1 = (-l) k -co 2k -A, 



(23) 



whence from (16) it is immediate that [letting 7 denote 7 2 ] 
exp(A-§) = exp 



o , /V 




-co 2 7-4), 

- 





cos(g)§)I, [l/oo]-5m(oo(|))-7 
-Q)-sin((a-$)I, cos((d-$)I 



(24) 



Principal Result 



THEOREM . Consider the semi-restricted Copernican-Newtonian (n + \)-body problem as formulated 
in (7a) and (9), with Initial Conditions (ICs) corresponding to Keplerian concentric circular 
orbits, i.e. in which 

/(0)-v k (0) = 0, (25a) 

so that 3 parameters suffice to specify the ICs. Let these parameters be defined by the triplet 
(p,9,co) given by polar coordinates in which p denotes the initial orbital radius, 9 the initial 
phase angle, and CO the initial angular velocity, assumed to conform to Kepler's Law (8), which 
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reduces the IC to a 2-parameter set (p,9). Specifically, assume (8) and let the ICs of (7a) & (9) 
be 

x\0) = p(cos(e\sin(e°)y , v k (0) = CO -J -x\0) = - ® p(cos(Q°),sin(Q°)y, (25b) 

k k k k 2 k k k k 

where ' denotes vector-matrix transposition. Suppose further that the ICs are so chosen that the 
frequencies are resonant as in (1); for later convenience, we may assume that 

1 < m, m < 5, (26) 

i j 

although the main result holds when these integers are arbitrarily large, which means that the 

initial frequencies (cfl.,Cfl.) are essentially arbitrary. Now, keeping the resonant frequencies 

(cfl.,co.) found at p = fixed, consider the variation of (p,9) as a function of p in such a way as 

tO preserve the periodicity [namely, "isoperiodic" continuation]: 

x\t + 7) = x\t), T = 2-71/co, co = co./m. = co./m., (27) 

while ignoring the question of periodicity of x ] , j = i + 1. A necessary and sufficient condition 
that for sufficiently small |i there exist ( p.(|l), 6(|l) ) preserving the periodicity (27) is that 
there exist a phase-shift § = (j)((i) and a distal multiplier (3 = (3(|l) such that 

Q m Ql) = G.Oi) + 0, (28) 

p. +1 0i) = P-p.(^i), (29) 
where, in the second post-Keplerian approximation ((|),(3) are given by 

<|> = m-(n/2) + ( m = 1, 2, 3, - ), (30a) 

+ - = 1/J 

whence, by Kepler's Law in the form 



p = P + - = lAl(3/2) 2/3 - 1 + - = 1.794980 + - , (30b) 



r = m/m = p 3/2 + - = 2.40 + (30c) 

i i+1 

the ONLY [!] possibilities for the resonance in the low-order case (26) are 

m:m = 2:1 or 5:2 or 3:1. (30d) 
i i+i v ' 

REMARK . Cancel the arbitrary supplementary assumption (26) which leads to (30d); then the 
initial generating orbits can be placed arbitrarily closely to ANY pair of concentric circles with 
arbitrary radii (because an arbitrary irrational P in (29), and so its (3/2) th power r as in (30c), 
can be approximated as closely as desired by a rational number m/m,^). In this way it can be seen 
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that the unique distal multiplier (3 given in (30b) as 1.80 at u: = is a constant which 
is completely independent of the assumed trial value of (3 at |i = ! 



PROOF . Let e denote the columns of 1^ = (e l ,e 2 ). Then at fl = the solution of the given 
ODE ICP is 

x\t) = p expiJ -0 )-e, = co -t + 9° = m o + 9°, o := a t. (31) 

k 2 k k kk k k 

We want this to be the linear part of the reformulation presented in the Lemma, so choose A as in 

(24), and then re-express the problem in the equivalent form (15), where the required matrix 

exponential is given by (24); here now z = Ot',v')' e E 4 and, correspondingly, the g(z) in (15) is 
given by g = ( 0', [/ m ]')', where 

/ m = . . x Vj - [* m /ll* m ll 3 ])+[(co k ) 2 " G-d + ^ k )/ii/ii 3 ]-/. (32) 

Thus, define the 'osculating' generating solution at \i = by 

x k,0 (f) = co5(co t)-x k (0) + [1/co ]-5m(co 0-v k (0), (33a) 

k k k 

v k '°(t) = cos((£> -t)-v k (0) - CO sin((£> t) x k (0). (33b) 

k k k 

Now the problem is rigorously equivalent to solving the integral equation 

t 

x\t) = x k '°(t) + f [1/co ]-sm(co [t - x])-g km (x k (x),x m (x)) dx; (34a) 

J k k 



although the formulation (24) gives a second equation 

t 

v k (t) = v k '°(t) + J cos(co k [? - x]yg km (x k (x),x m (x)) dx, (34b) 



the second is a consequence of differentiating the first with respect to time t, and therefore is of 
no further consequence, after noting that in (33a) the initial velocity v k (0) = C0 k -7 2 ^ k (0) is 
already defined in terms of ^ k (0) by hypothesis. 

It is a well-known consequence of basic ODE theory that the (necessarily unique) solutions of 
(7a) and (9) for |i = 0, namely (33), and the solutions for sufficiently small |i > remain 
arbitrarily close together for any finite time, specifically here for < t < T = 2-7i/co, provided 
only that |i be sufficiently small. Therefore for sufficiently small |i one knows a priori that the 
solutions remain within planar concentric annuli surrounding the circular orbits (33), and therefore 
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do not approach each other during the time of interest. Accordingly one may find a global Lipschitz 
constant for such relevant domains of the Cartesian products of E 4 (in which the Jacobian matrix of 
the right-hand side of (11) is continuous, and so has a bounded norm K in the relevant z-domain). 
Therefore the Lemma is applicable, for sufficiently small [i, and so the solution of (34) exists and 

can be constructed by successive approximations (with j = 1, 2, 3, ■■■) : 

t 

x k '\t) = x k '°(t) + f [l/co ]sm(co [t - T])-g% Xj ~\T),x m ' H (T)) dx, (35) 

J k k 



taking as the initial iterate the decoupled Keplerian solution (33), for < t < T. 

Next, look at the solution of (34) only at every revolution of duration T = 2-7i/co, wherein the 
commensurability of co and co. +i ensures that there is a common period to the initial iterates of the 
two adjacent orbits. This give the famous Poincare map: 

271 

x\T) = x k (0) - Lt-(GLi°/[m co 2 ])- \sin(m a)-h km (x\a),x m (a)) da, (36a) 

m k J k 



h km := [[( X m - x k )/\\x m - x k ll 3 ] - [x m /llx m ll 3 ]] + [(a/ - Gp (l + u^/ii/h V 
We shall prove that this map can be expressed in the form 



(36b) 



x\T) = x\0) - |l-(G|l°/[my])/ m (p 0°p ,9°;|l), (37a) 

m k k k m m 



271 

/ m (p 6°p ,6°;|i) := f™(mo)AW>)) da, (37b) 

k k m m J k 



where for certain specific values of (P k >6°) the Jacobian determinant of /""with respect to (P k >9°) 
is positive at fl = 0, and where = for these same specific values, at |i = 0; then, by the 
Implicit Function Theorem there exist functions (p k (|i),9 k (|l)) for sufficiently small values of |i 
which satisfy /™ = and so which provide the desired periodicity of x k {T). 

It will become evident that the conditions to be derived are also necessary, because if they 
are not satisfied, provides a resonant forcing term which drives x k (T) ever farther from x k (0), 
and will do so without bound (or until the resonant 'pumping' drives the map's next iterate out of 
the map's domain of definition). 

To anticipate the results of a rather lengthy and arduous calculation, it will be proved below 
that y 1 ™ has the form 

/ m = r| exp(J O )vc m (0) - r| exp(J O ) /(0), (38) 

J 'km,m ^ V 2 km,rn V ' 'km,k ^ V 2 km.k 7 V " V ' 

where the scalar functions T| k (|i) and <I> k (|i) have the properties that 
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r| (0) = 0, r| (0) = 1/2, O (0) = 7i/2, (39) 

'km,m V ' ' 'km,k V ' ' km,k V ' y ' 

so that at fl = the Jacobian matrix H of with respect to x k (0) is given by 

H = -(l/2yexp(J 2 [n/2]) = -(1/2)-/ <fe/(fl) = 1/4 > 0, (40) 

and, as claimed, the Leray-Schauder Index [29] of the chosen generating solution is unity. 

It is evident by inspection that upon the second iteration the third term (36b) drops out by 
the choice of a Keplerian generating solution, so that for consideration of the second iterate of 
the successive approximations we may simplify (36b) to 

; km r/ m kw.. m k.^n r m /M m.,3-. /a-i\ 

h := [(x - x )/\\x - x II ] - [x / II jc II J. (41) 
It is also immediate that the second term in (41) may be omitted, because by (31) 

x m /\\x m \\ 3 = (l/p 2 )-exp(J m c) exp(J Q°) e l , sin(m c) exp(J m c) do = (42) 

m 2m 2m J k 2m 



because of the hypothesis that m £ m . Thus we have simplified the calculation to that of 

km 

evaluation of 

2TC 

/ m = \sin(m o)- , P(o)-(/(o) - /(c)) do, x\o) = p exp(J [m o + 9 ]) *? 1 , (43) 

J k k 2 k k 



:= l/\\x m (o) - x\o)\\ 3 . (44) 

Now by the chief property of the scalar product (x,y), namely that if M is an arbitrary matrix 
(Mx,y) = (x,M'y), and the fact that / ' = - / , it is easy to calculate that 

II x m (o) - x\o) II 2 = lljc m ll 2 + HxV - 2-p ■ p k (e l ,exp(jpe l ), (45a) 
t = (m - m >o + 9° - 9°, (45b) 

^ V k m k m' ' 

V(o) = i|/(o)/[ (p k ) 2 + (p m ) 2 ] 3/2 , (45c) 
V (o) = l/[ 1 - (zJ 2 cos(k) f'\ (ej 2 := 2-p^pJf (p/ + (pj 2 ], (45d) 
\|/(o) = 1 + (3/2)-(e ) 2 co5(c;) + (15/8)-(e ) 4 co/(c;) + ■■■ , (45e) 

km km 

where the binomial series in (45e) always converges because (e ) < 1 whenever p and p are 

km k m 

2 

distinct (which is a consequence of the fact that then < (p k - p ) and of obvious 



manipulations of the expansion of the latter). 
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Now, remembering that \|/ depends upon both k and m, define 

271 

A := (1/[2ji])- \ sinim o)-cos(mo)\\f(o) do, 

kmj J k j 





271 



5 := (l/[2n]Y sm(m o)-sinffflo)f(o) Jo, 

kmj J k j 



kmj 



{ (A ) 2 + (5 ) 2 } 1/2 , 

1 V krn V km 7 J ' 



(46a) 
(46b) 

(46c) 



O := -Arctan{ B IA }, 

kmj kmj kmj 



and note that (43)-(44) may be simplified by use of the novel identity 

271 

(l/[27i]> \sin(m 6)-\u(6)-exp(J [mc + 9°]) do = ri exp(J [O +9°]). 

J k 2 j j kmj 2 mkj j 



(46d) 



(47) 



It is the radical simplification provided by the apparently hitherto unnoticed identity (47) which 
appears to be the chief innovation in the present work. For we may now write that 



fn 



d^/[ (p/ + (pj 2 f 2 , 



d km , 



p r| exp(JJ® + 9°]) - p Tl exp(JJ® + 9°]) 

^m 'kmm ^ v 2 L kmm m J/ ^k 'kmk ^ v 2 L kmk k J/ 



■e = 



(48a) 



= p Tl 

m kmm 



cos(0 +9°) 

kmm m 

sin(® +9°) 

kmm m 



- P Tl 

^k 'kmk 



co5(0 + 9°) 

kmk k 

sin(Q + 9°) 

v kmk r 



(48b) 



From mere inspection of (48) it is now evident that the necessary and sufficient conditions for 
y 1 ™ to vanish are that there be an orbital resonance defined by 



P m = P-Pv 



B = Ti m 

kmk kmm 



and a corresponding enabling phase-shift 



9° = 9° + d) + M-n, $ = O - O , ( M = 1, 2, 3, - ). 

m k mk mk kmk kmm 



(49a) 



(49b) 



This completes the easy part of the present derivation. 

Now begins the hard work of evaluation of B and (|). Part of this is easy, because by inspection 
we need only to evaluate the lowest-order terms in 



A = + - 

kmk 



B = (1/2) + 

kmk v ' 



(50) 
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because in this case it is adequate to use \|/ = 1 + ■■■ because of the well known orthogonality 
properties of sines and cosines. In contrast, the lowest order terms in A and B vanish 

kmm kmm 

identically, and we must go to the second-order terms in \|/ in order to get meaningful results. Thus, 
to the second order in the series (45e) 

(A , B ) = (3/2)(£ ) 2 ( a , (3 ) + ■■■, (51a) 

v kmm kmm' v ' v km v kmm 'kmm ' ' 

271 

a := (l/[2jt]Y ( sinim c)-cos(m o)-cos([m - m ] o + 9° - 9°) do, (51b) 

kmm J k m km km 



271 

(3 := (l/[27i]> \sin(m c)-sin(m o)-cos([m - m ] o + 9° - 9°) do, (51c) 

kmm J k m km km 



which is where conceptualizing ends and labor begins. My advice to the reader is to replace each 
sine and cosine by the sum or difference of two complex exponentials, as in de Moivre's Theorem, and 
then multiply out the resulting 6 products as complex numbers. From this lengthy exercise in 
elementary complex algebra there results the [ultimately real] identities: 

a = -(l/4)-sm(9° - 9°), (3 = (l/4)-cos(9° - 9°), (51d) 

kmm V ' V k rn K kmm V ' V k rn V ' 

2 2 1/2 

from which we obtain the welcome simplification that [ (a ) + (B ) 1 = 1/4. Therefore, 

kmm kmm 

finally, by (46c) 

r| = (3/2>(e ) 2 (l/4) + -, r| = (1/2) + (51c) 

'kmm V ' V kxa V ' ' 'kmk \ ) ■> V J 

so that, by (49a), to lowest order in |i, 



1/(3 = p /p = n lx\ = (3/2) p p /[ (p ) 2 + (p ) 2 ] 3/2 

K K k K m 'kmm 'kmk V ' K k K m L VK k 7 yV m J 

= (3/2)-(l/(3)/[ 1 + (l/(3) 2 f 2 , (51f) 

which requires for self-consistency that, reminiscent of the 1766/1772 Titius/Bode 'Law', 

1 + (l/(3) 2 = (3/2) 2/3 , (51 g ) 

i.e. that (3 have the unique particular value claimed in (30b)! Note that to lowest order this (3 is 
a universal constant, independent of the gravitational constant G or the masses of the planets! 
The proof of (30a) is analogous but simpler, noting that from (46d) and (5 Id) 

(j) = O - O = -(9° - 9°), (51h) 

mk kmk kmm m k 

so that, bringing (9° - 9°) to the left-hand side of (49b) and dividing by 2 we obtain the claimed 

m k 

result (30a). Equation (38) is a trivial consequence of (48b). Also evaluation of the Jacobian with respect to (p, 9) 

instead of with respect to the components of makes x(0) makes no difference to the claim that the Leray-Schauder 
Index of the generating solution isolated by (51g) is unity. This completes the proof. 



Now that the importance of the distal multiplier (3 for the distance beyond an inner planet for 
an outer planet in terms of the ratio of orbital radii has been derived in full rigor for the 
coplanar, concentric 3-body problem, it seems permissible to use a cruder physical model to consider 
the case of 4 or more bodies. For simplicity, keep the first body anchored at the origin and let 
the bodies have mean orbital radii p ; (z = 1, 2, ■■■ n), where p L = by definition, and 

= Pi < P 2 < •" < Pi < p i+1 < - < p n , (55) 

and where each body (initially assumed decoupled from mutual interactions) is started with Keplerian 
circular velocity v ; = Vg-M/ J~p i where M = m x denotes the large central mass and where the smaller 
masses are denoted by m i = uyM in terms of ratios « 1 (i = 2, 3, ■■■ , n). Let r o > denote 
initial total angular momentum, and E = 1 2^ | = -Eq > denote the absolute value of the initial 
total energy. Also, define y = T^Ig-M 3 ' 2 and e = 2-EJ^G-M 2 , and it is easy to verify that 
(neglecting terms quadratic or higher in the jO-j) Conservation of Angular Momentum and Conservation 
of Energy (combined with the Virial Theorem) are given by y = y and e = e where 



£ = 



v m/pi 



i=2 



L 

i=2 



(56) 



THEOREM . If the planetary orbits satisfy a Titius-Bode Law of the form 

Pi = |r 2 -p 2 , (i = 2, 3, 4, - , n), (57) 

then the distal ratio (3 = z 2 is the square of the unique positive root z > of a polynomial of 
degree 4-n - 8 in z whose coefficients are functions only of the constants y and e and the mass- 
ratios (ij and which can be derived by elimination of p 2 between the expressions for y and e in (56). 

PROOF . Insert (57) into (56) and compute y 2 £. I shall publish a general algorithm defining 
all of the coefficients explicitly elsewhere; however, it can be recovered easily by the interested 
reader after following the next example in the case n = 3, wherein the algebra is less difficult. □ 

REMARK . To apply this result to the present solar system, simply replace the p { in (56) by the 
actual mean distances R { from observational astronomy and use the actual mass-ratios ja i from 
astrophysics in order to evaluate y and e numerically. After finding (3, the eliminated p 2 can be 
recovered by solving either y = y or 8 = e following which the remaining p { are given by (57). 
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In the case n = 3, the polynomial F(z) whose roots are sought has the form 

F(z) = z 4 + a 3 -z 3 + cl 2 -z 2 + a^z 3 + oc , (58) 

which, as in the proof of the preceding theorem, can be obtained simply by comparing coefficients in 
the multiplied-out version of 

(|i f z 2 + u\ 2 Hm + [i-izf - az 2 = 0, a = e -(y ) 2 . (59) 

The result is 

a 3 = 2-(fi 1 /u, 2 ), a 2 = (n^) 2 + (M-2/M-1) - a/iii^) 2 ], ^ = 2, a = (i/^. (60) 

Using data for Jupiter and Saturn, the reader may verify that F = has only one positive root, 
whosee square yields (3 = 1.833. With more effort, analogous results can be obtained for arbitrary 
n (see Table 1 ). Note that the value ofb had converged to 1.795, in amazing agreement with the value obtained in 
(30b). 



m { 




Pi 


BODY 


n 


P 


Pi 


1,000,000 






Sun 


1 






0.16601 


0.3871 


0.2818 


+Mercury 


2 






2.447841 


0.7233 


0.50608 


+Venus 


3 


1.86 


1.868 


3.003469 


1.000 


0.9088 


+Earth 


4 


1.46 


1.382 


0.322714 


1.5237 


1.632 


+Mars 


5 


1.47 


1.5237 


A 








6 






954.60258 


5.2030 


5.26 


+Jupiter 


7 


1.68 


1.927 


285.80765 


9.5281 


9.45 


+Saturn 


8 


1.766 


1.831 


43.549846 


19.1829 


16.97 


+Uranus 


9 


1.800 


2.013 


51.67161 


30.0796 


30.48 


+Neptune 


10 


1.795 


1.568 


0.007541 


49.0250 


54.70 


+Pluto 


11 


1.795 


1.6298 



Table 1 
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Conclusions 

In this paper, we showed that in a planetary or satellite system, resonant perturbations will cause 
dynamical evolution of the orbital radii except under highly specific mutual relationships. Simplifying 
the problem by reformulating it as a hierarchical sequence of [unrestricted] 3-body problems, in which 
gravitational interactions are ignored except between the central body and the body of interest and the 
next outwardly orbitally adjacent body, it is proved that the resonant perturbations from the outer body 
will destabilize the inner body [& conversely] unless its mean orbital radius is a unique and specific 
multiple (3 of that of the inner body . In this way a sequence of concentric orbits can each stabilize the 
next adacent inner orbit, until only the outermost orbit remains; but it is already tied to the collection of 
inner orbits by conservation of total angular momentum and so the entire configuration becomes stabilized. 
Let fj, = M Max /M denote the ratio of the mass of the largest small body M M ax to that, M,of the large 
central body; in our Solar System, \i is less than 1CT 3 . Expanding (3 in a power series in /i, the lowest-order 
terms for distal multiplier (3 and phase shift are found to start with the following universal constants 
[for m = 1,2,3] 

(3 = (3 + = l/^/(3/2)(2/3) _ i = 1.794980 + , = m(ir/2) + .. 

which agrees with the observed distal ratios between Jupiter and the Asteroids, and between Saturn & 
Jupiter, and Uranus & Saturn. 
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